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We study asymptotically constrained systems for numerical integration of the Einstein equations, 
which are intended to be robust against perturbative errors for the free evolution of the initial data. 
First, we examine the previously proposed "A-system", which introduces artificial flows to constraint 
surfaces based on the symmetric hyperbolic formulation. We show that this system works as expected 
for the wave propagation problem in the Maxwell system and in general relativity using Ashtekar's 
connection formulation. Second, we propose a new mechanism to control the stability, which we 
call the "adjusted system". This is simply obtained by adding constraint terms in the dynamical 
equations and adjusting its multipliers. We explain why a particular choice of multiplier reduces 
the numerical errors from non-positive or pure-imaginary eigenvalues of the adjusted constraint 
propagation equations. This "adjusted system" is also tested in the Maxwell system and in the 
Ashtekar's system. This mechanism affects more than the system's symmetric hyperbolicity. 

PACS numbers: 04.20.Cv 04.25.-g 04.25.Dm 



I. INTRODUCTION 



Numerical relativity, an approach to solve the Einstein equations numerically, is supposed to be the only way to 
study highly non-linear gravitational phenomena. Although the attempt has already decades of history, we still do 
not have a definite recipe for integrating the Einstein equations that will give us us accurate and long-term stable 
time evolutions. Here and hereafter, we mean "stable evolution" that the system keeps the violation of the constraints 
within a suitable small value in its free numerical evolution. 

As the authors discussed in our preceding paper (Paper I) Q, one direction for obtaining a more stable system 
is to apply a set of dynamical equations which have manifest hyperbolic form (or first-order form). The standard 
Arnowitt-Deser-Misner (ADM) formulation does not have this feature, but there are many alternative proposals for 
constructing a hyperbolic set of equations ( |^-^|, see also references in However, we showed in Paper I that a 
symmetric hyperbolic form (mathematically, the "ultimate" level of hyperbolicity) does not necessary give the best 
performance for stable numerical evolution compared with weakly and strongly hyperbolic systems. This experiment 
was performed using Ashtekar's connection variables 0, since this formulation enables us to compare three levels of 
hyperbolic formulations keeping the same fundamental dynamical variables. 

In this article, we discuss different (but somewhat related) approaches to obtaining stable evolution of the Einstein 
equations. The idea is to construct a system robust against the perturbative error produced during numerical time 
integration. We discuss the following two systems. 

The first one is the so-called "A-system" , which was proposed originally by Brodbeck, Frittelli, Hubner and Reula 
(BFHR) ||. The idea of this approach is to introduce additional variables, A, which indicates the violation of 
the constraints, and to construct a symmetric hyperbolic system for both the original variables and A together 
with imposing dissipative dynamical equations for As. BFHR constructed their A-system based on Frittelli-Reula's 
symmetric hyperbolic formulation of the Einstein equations [|J, and we |)| have also presented a similar system for 
Ashtekar's connection formulation [Q] based on its symmetric hyperbolic expression [ [lC| jTT[ ] . In §0, we review this 
system and present numerical examples which show this system behaves as expected. 

The second one has the same motivation but turns to be more practical, which we call "adjusted-system" . The 
essential procedure is to to add constraint terms to the right-hand-side of the dynamical equations with multipliers, 
and to choose the multipliers so as to decrease the violation of the constraint equations. This second step will be 
explained by obtaining non-positive or pure-imaginary eigenvalues of the adjusted constraint propagation equations. 
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We remark that adjusting the dynamical equation using the constraints is not a new idea. This can be seen for 
example in a remedial ADM system by Detweiler , in a conformally decoupled trace-free re-formulation of ADM 
by Nakamura et al jL3|, and also in constructing hyperbolic formulations We also remark that this eigenvalue 



criterion is also the core part of the theoretical support of the above A-system. In §111, we describe this approach and 



present numerical examples again in the Maxwell system and in the Ashtekar's system. 

This "adjusted-system" does not change the number of dynamical variables, and does not require hyperbolicity in 
the original set of equations. Therefore we think our results promote further applications in numerical relativity. 

We will not repeat our explanation of Ashtekar's connection formulation in our notation, nor of our detailed 
numerical procedures, since they are described in our Paper I §. 



II. ASYMPTOTICALLY CONSTRAINED SYSTEM 1: A-SYSTEM 



We begin by reviewing the fundamental procedures of the "A-system" proposed by Brodbeck, Frittelli, Hubner and 
Reula (BFHR) ||. We, then, demonstrate how this system works in Maxwell's equations, and Ashtekar's connection 
formulation of the Einstein equations in the following subsections. 



A. The "A system" 

The actual procedures for constructing a A system are followings. 

(1) Prepare a symmetric hyperbolic evolution system which describe the problem; say 

d t v? = A* 6 8 i v, s + B r ', (2.1) 

where u 1 (7 = 1,..., N) is a set of dynamical variables, A(u(x 1 )) forms a symmetric matrix (Hermitian matrix 
when u is complex variables) and B(u(x 1 )) is a vector, where A and B do not include any further spatial 
derivatives in these components. The system may have constraint equations, which should be the first class. 
Ideally, we expect that the evolution equation of the set of constraints C p (p = 1, . . . , M), which hereafter we 
denote constraint propagation equation, forms a first order hyperbolic system (cf. ]l4[]), say 

d t C p = D^adiC + E^C 7 , (2.2) 



(where D, E are the same with A, B above) but this hyperbolicity may not be necessary. 

(2) Introduce A p as a measure of violation of the constraint equation, C p ~ 0. (« denotes "weakly equal".) Here 
C p is a given function of u and is assumed to be linear in its first-order space derivatives. We impose that A 
obeys a dissipative equation of motion 

d t X p — ct( p )C p — /?( P )A P (we do not sum over p and (p) on right hand side) (2.3) 



with the initial data A p = 0, and by setting a ^ 0, (3 > 0. We remark that A p remains zero during the time 
evolution if there is no violation of the constraints. 

(3) Take a set (u, A) of dynamical variables, and modify the evolution equations so as to form a symmetric hyperbolic 
system. That is, the set of equations 



(= means that we have extracted only the term which appears in the principal part of the system) can be 
modified as 

*0?W£: VMS). m 



where the additional terms will not disturb the hyperbolicity of equations of u 1 , rather they make the whole 
system symmetric hyperbolic, which guarantees the well-posedness of the system. 
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Therefore the derived system, (2.5), should have unique solution. If a perturbative violation of constraints, X p ^= 0, 
occurs during the evolution, by choosing appropriate as and /3s in (|]^), As can be made decaying to zero, which means 
the total system evolves into the constraint surface asymptotically. We note that this procedure requires that the 
original system u forms a symmetric hyperbolic system, so that applications to the Einstein equations are somewhat 
restricted. BFHR M constructed this A-system using a Frittelli-Rcula's formulation We || also applied this 
system to the symmetric hyperbolic version of Ashtekar's formulation [jlp|l . 

We next review a brief proof why the system ( |2.5| ) ensures that the evolution is constrained asymptotically. We 
first remark again that we only consider perturbative violations of constraints in our evolving system. The steps are 
following. 

(a) Since we modify the equations for it 7 , the propagation equation of the constraints are also modified; write them 
schematically as 

d t C p = D'P^C 7 + E p a C a + G 1311 ^^ + H lp a dX + I p a X a . (2.6) 

(b) In order to see the asymptotic behaviors of (A p , C p ), we write them using their Fourier components so that their 
evolution equations take an homogenous form. That is, we transform (A P ,C P ) to (X P ,C P ) as 



X(x,t) p = X(k,t) p exp(ik- x)d i k, C(x,t) p = C(k,t) p exp(ik ■ x)d 6 k. (2.7) 



Then we see the evolution equations (2J3) and (2_i) become 



a{p)K v ? ) =■ p n° i (28) 



(c) If all eigenvalues of this coefficient matrix P have negative real part, a pair (A, C) evolves as exp(— At) asymp- 
totically where —A is the diagonalized matrix of P, which indicates that the original variables (A, C) evolves 
similarly. It would be best if we could determine the a and f3 in such a way in general, but it is not possible. 
Therefore we extract the principal order of P and examine the condition for a and /3 so that P only has negative 
(real) eigenvalues. We remark again that this procedure is justified when we only consider a perturbative error 
from the constraint surface. 



B. Example 1: Maxwell equations 



As a first example, we present the Maxwell equations in a form of A-system. The Maxwell equations form linear and 
symmetric hyperbolic dynamical equations, together with two constraint equations, which might be the best system 
to start with. 



1. X-system 

The Maxwell equations for an electric field E l and a magnetic field B l in the vacuum consist of two constraint 
equations, 

C E := « 0, (2.9) 
C B := d,B l « 0, (2.10) 

and a set of dynamical equations, 




(2.11) 



which satisfies symmetric hyperbolicity. The constraint evolutions become c^Ce = and dtCs — 0, which indicate 
(trivial) symmetric hyperbolicity. According to the above procedure, we introduce As which obey 
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2. Analysis of eigenvalues 

Now the evolution equations for the constraints Ce and Cb become 

d t C E = ai(AA E ), d t C B = a 2 (AX B ) 



(2.16) 
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where A = did 1 . We take the Fourier integrals for constraints Cs [(2.16)] and As [( |2.12| ), (2.13)], in the form of (|2.7|). 
to obtain 



(2.17) 



where k 2 — kjk\ W e find the matrix is constant . Note that this is exact expression. Since the eigenvalues are 
(—pi ± v/3f — 4a 2 fc 2 )/2 and (— /3 2 ± \//?| — 4c*2fc 2 )/2, the negative eigenvalue requirement becomes ai,a 2 ^ and 

01,02>O. 



3. Numerical demonstration 



We present a numerical demonstration of the above Maxwell "A-system" . We prepare a code which produces elec- 
tromagnetic propagation in xy-plane, and monitor the violation of the constraint during time integration. Specifically 
we prepare the initial data with a Gaussian packet at the origin, 



E l (x,y,z) = (-Aye- B{x2+v2) ,Axe- B ^ 2+ y 2 \0), 
B l (x,y,z) = (0,0,0), 



(2.18) 
(2.19) 



where A and B are constants, and let it propagate freely, under the periodic boundary condition. 

The code itself is quite stable for this problem. In Fig.|l|, we plot L2 norm of the error (Ce over the whole grid) 
as a function of time. The solid line (constant) in Fig.|l| (a) is of the original Maxwell equation. If we introduce As, 
then we see the error will be reduced by a particular choice of a and j3. Fig|l] (a) is for changing a with (3 — 2.0, 
while Figjl] (b) is for changing (3 with a = 0.5. Here, we simply use a := a% = a 2 and (3 := j3\ = (3 2 - We see better 
performance for (3 > [Fig.|l| (b)], which is the case of negative eigenvalues of the constraint propagation equation. We 
also see the system will diverge for large a [Fig.|l| (b)] . The upper bound of a can be explained by the violation of the 
Courant-F riedri ch-Lewy (CFL) condition, where the characteristic speed comes from the flux term of the dynamical 
equations ( 2.15 ). 
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FIG. 1. Demonstration of the A-system in the Maxwell equation. Fig. (a) is constraint violation (L2 norm of Ce) versus 
time with constant /3(= 2.0) but changing a. Here a = means no A-system. Fig.(b) is the same plot with constant a(= 0.5) 
but changing (3. We see bet ter p erformance for f3 > 0, which is the case of negative eigenvalues of the constraint propagation 
equation. The constants in (2.18) were chosen as A — 200 and B — 1. 



C. Example 2: Einstein equations (Ashtekar equations) 



The second demonstration is of the vacuum Einstein equations in Ashtekar's connection formalism 0. 

Before going through the A-system, we will briefly outline the equations. The fundamental Ashtekar variables are 
the densitized inverse triad, E l a , and the SO(3,C) self-dual connection, Af, where the indices indicate the 

3-spacetime, and a, 6, ••• is for SO(3) space. The total four-dimensional spacetime is described together with the 
gauge variables N, N l , Aq, which we call the densitized lapse function, shift vector and the triad lapse function. Since 
the Hilbert action takes the form 



d 4 x[(dtA2)K + NC H + fPCm + AqCgo], 



(2.20) 



the system has three constraint equations, Cr ~ Cmi ~ Cc a ~ 0, which are called the Hamiltonian, momentum, and 
Gauss constraint equation, respectively. They are written as 



where 
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equation for (E^,Af) constitutes a weakly hyperbolic form, 



(2.21) 
(2.22) 
(2.23) 

diE 3 a — ie a b c A^El. The original dynamical 
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d t A1 = ~ie ab c NE° b F?j + N 3 F£ + V,A a a 



(2.24) 
(2.25) 



where X>,AT := d s Xg -Je^A b j X J c i , for X l J + Xg = 0. It is also possible to express ( |2.24| ) and ( |2.25| ) to reveal 



symmetric hyperbolicity |lO|]ll[ | . For more detailed definitions and our notation, please see Appendix A of our Paper 

10- 
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1. X-system for controlling constraint violations 



Here, we only consider the A-system which controls the violation of the constraint equations. In ||, we have also 
discussed an advanced version of the A-system which controls the violations of the reality condition. 
We introduce new variables (A, Ai, A a ), obeying the dissipative evolution equations, 

d t X = a 1 C H - ft A, (2.26) 
d t \ i =a 2 C m -fa\ ii (2.27) 
d t \ a = a 3 C Ga ~ (3 3 \ a , (2.28) 

where Oj 7^ (possibly complex) and fit > (real) are constants. 



If we take y a :— (E^,Af, A, A,, A a ) as a set of dynamical variables, then the principal part of ( |2.26| )-( |2.28| ) can be 
written as 



dt\i^a 2 [-e6 l i Ei + e5iE l b ](d l A<}), 
d t X a = a 3 diE l a . 



(2.29) 
(2.30) 
(2.31) 



The characteristic matrix of the system u a is not Hermitian. However, if we modify the right-hand-side of the evo- 
lution equation of (E a , A%), then the set becomes a symmetric hyperbolic system. This is done by adding CI37 (9;A ) 
to the equation of d t E a , and by adding ia 1 e a c d E^E l d (diX) + a 2 {-e^ lm Ef + e5™E la )(di\ m ) to the equation of d t Af. 
The final principal part, then, is written as 
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(2.32) 



(2.33) 



(2.34) 



Clearly, the solution (E a , Af, A, A,-, A a ) = (E a , Af, 0, 0, 0) represents the original solution of the Ashtekar system. 



2. Analysis of eigenvalues 



After linearizing and taking the Fourier transformation (2.7), the propagation equation of the constraints 
(C H ,CMi,C Ga ) and (A, A;,A ) are written as, 
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(2.35) 



In order to link the discussion with our later demonstration in the plane symmetric spacetime, we here consider 
only the F ourie r component of fej = (1, 0, 0) for simplicity. The eigenvalues, Ei (i = 1, • • • , 14), of the characteristic 
matrix of (2.35) can be written explicitly as 



G 



(E u - ■ ■ , E w ) = -(1/2)03 ± (1/2)^| - 4|a 3 | 2 , 

-(l/2)(i + ft) ± (1/2)^-1 -4|a 3 | 2 -2*ft + ft 2 , 

-(l/2)(-i + ft) ± (1/2)^-1 -4|a 3 | 2 -2ift. + /?|, 

-(l/2)(t + ft) ± (1/2)^-1 -4|a 2 | 2 -2ift + /?f, 

-(l/2)(-i + ft) ± (1/2)^-1 - 4|« 2 p - 2i/3 2 + /5| 

and as solutions (En, ■ • • , -E14) of the quartic equation 

+ (ft + ft)x 3 + (2| ai | 2 + 2\a 2 \ 2 + 1 + ftft)x 2 + (2|a 2 | 2 ft + ft + ft + 2| ai | 2 ft)x + (ftft + 4|a 1 | 2 |a 2 | 2 ) = 0, 

(2.36) 

where |aq| 2 = a^a^. We omit the explicit expressions of En, ■ ■ ■ , E14 in order to save space. 
A possible set of conditions on a p , ft, (p = 1, 2, 3) for die(Ei) < are 

a p ^ and ft > 0. (2.37) 

This is true (necessary and sufficient) for Ei, ■ ■ ■ , E±q, and also plausible for En, ■ ■ ■ , E\± as far as our numerical 
evaluation tells (see Fig.|f). 
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FIG. 2. Example of eigenvalues of the system ( |2.35| ). We plot the eigenvalue which has the maximum real part between 
four of them, for the case of fixing ai — Q2 = 1 and changing ft and ft. We see our desired condition, "all negative real 
eigenvalues", is available when the combinations produce the solid lines. That is, when both /3 take the large positive values. 



3. Numerical demonstration 



In this subsection, we demonstrate that the A-system for the Ashtekar equation actually works as expected. 

The model we present here is gravitational wave propagation in a planar spacetime under periodic boundary 
condition. We perform a full numerical simulation using Ashtekar's variables. We prepare two +-mode strong pulse 
waves initially by solving the ADM Hamiltonian constraint equation, using York-O'Murchadha's conformal approach. 
Then we transform the initial Cauchy data (3-metric and extrinsic curvature) into the connection variables, (E^,Af), 
and evolve them using the dynamical equations. For the presentation in this article, we apply the geodesic slicing 
condition (ADM lapse N = 1, with zero shift and zero triad lapse). We have used both the Brailovskaya integration 
scheme, which is a second order predictor-corrector method, and the so-called iterative Crank-Nicholson integration 
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scheme for numerical time evolutions. The details of the numerical method are described in the Paper I [|J, where we 
also described how our code shows second order convergence behaviour. 

In order to show the expected "stabilization behaviour" clearly, we artificially add an error in the middle of the 
time evolution. More specifically, we set our initial guess 3-metric as 





sym. 1- K{e~( x - L f +e- { - x+L ?) : 



in the periodically bounded region x = [—5, +5], and added an artificial inconsistent rescaling once at time t — 6 for 
the Ay component as A y — ► Ay(l + error). Here K and L are constants and we set K = 0.3 and L = 2.5 for the plots. 




FIG. 3. Demonstration of the A-system in the Ashtekar equation. We plot the violation of the constraint (L2 norm of 
the Hamiltonian constraint equation, Ch ) for the cases of plane wave propagation under the periodic boundary. To see the 
effect more clearly, we added artificial error at t = 6. Fig. (a) shows how the system goes bad depending on the amplitude 
of artificial error. The error was of the form Ay — ■> Ay (I + error). All the lines are of the evolution of Ashtekar's original 
equation (no A-system). Fig. (b) shows the effect of A-system. All the lines are 20% error amplitude, but shows the difference 
of evolution equations. The solid line is for Ashtekar's original equation (the same as in Fig. (a)), the dotted line is for the 
strongly hyperbolic Ashtekar's equation. Other lines are of A-systems, which produces better performance than that of the 
strongly hyperbolic system. 

Fig.^ (a) shows how the violation of the Hamiltonian constraint equation, Ch, become worse depending on the term 
error. The oscillation of the L2 norm Ch in the figure due to the pulse waves collide periodically in the numerical 
region. We, then, fix the error term as a 20% spike, and try to evolve the same data in different equations of motion, 
i.e., the original Ashtekar's equation [solid line in FigJ| (b)], strongly hyperbolic version of Ashtekar's equation (dotted 
line) and the above A-system equation (other lines) with different (3s but the same a. As we expected, all the A-system 
cases result in reducing the Hamiltonian constraint errors. 



D. Remarks for the A-system 



In the previous subsections, we showed that A-system works as we expected. The system evolves into a constraint 
surface asymptotically even if we added an error artificially. However, the A-system can not be introduced generally, 
because (i) the construction of the A-system requires that the original dynamical equation be in symmetric hyperbolic 
form, which is quite restrictive for the Einstein equations, (ii) the system requires many additional variables and we 
also need to evaluate all the constraint equations at every time steps, which is hard task in computation. Moreover, 
it is not clear that the A-system can control constraint equations which do not have any spatial differential terms. 
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(e.g., the primary metric reality condition in the Ashtekar formulation.) |^| 

We, next, propose an alternative system which also enable us to control the violation of constraint equations, but 
is robust for the above points. 

III. ASYMPTOTICALLY CONSTRAINED SYSTEM 2: ADJUSTED SYSTEM 

We here propose another approach for obtaining stable evolutions, which we name the "adjusted-system". The 
essential procedure is to add constraint terms to the right-hand-side of the dynamical equations with multipliers, and 
to choose the multipliers so as the adjusted equations decrease the violation of constraints during time evolution. This 
system has several advantages than the previous A-system. 

A. "Adjusted system" 

The actual procedure for constructing an adjusted system is as follows. 

(1) Prepare a set of evolution equations for dynamical variables and the first class constraints which describe the 
problem. It is not required that the system is in the first order form nor hyperbolic form. However here we 



start from the same form with ( |2.1| ) and (2.2). We repeat them as 

d tU ~< = A l ~< s d lU s + B-< , (3.1) 
8 t C p = D'p^C 7 + E p a C a , (3.2) 

where A(u(x 1 )) is not required to form a symmetric or Hermitian matrix. 



(2) Add the constraint terms, C p , (and/or their derivatives) to the dynamical equation (3.1) with multipliers k, 

dtu 1 = A^sdiU 6 + 5 7 + K JC + KpdiC. (3.3) 

We call the added terms, kJC p and/or K^diC 9 , "adjusted terms", and leave nj and kJ 1 unsp ecified for the 
moment. Because of these adjusted terms, the original constraint propagation equations, fl3.2| ), must also be 
adjusted: 

d t C = D ip a d t C a + E p a C° + F^ p <y d i d 3 C' J + G lp a diC a + H" ' a C° . (3.4) 
The last three terms are due to the adjusted terms. 



(3) Specify the multipliers k, b y ev aluating the eige nvalu es that appear in the RHS of (3.4). Practically, by taking 



the Fourier transformation (2.7), we can reduce (3.4) to homogeneous form, 

d t C p = (ihD ip a + E P (J - hkjF^^ + ihG lp a + H p a )&. (3.5) 

We, then, take a linearization against a certain background spacetime, 

dt {1) C P = (iki (0) L> ip CT + - kikj ^F ijp a + ih { ^G lp a + {0) H p a ) ^C a . (3.6) 



(here (n) indicates the order in linearization) and evaluate the eigenvalues of the coefficient matrix in (3.G). 
For this process, we propose two guidelines. 

(a) The first one is to obtain negative real-part of the eigenvalues. This is from the same principle in A-system 
when we specified a and (3, in order to force the system approach the constraint surface asymptotically. 
Provided that we obtain k which produce all the negative-real-part eigenvalues, the Fourier component C 
decays to zero in time evolution, and the original constraint term C also. 



1 This statement is not inconsistent with our previous work ||, in which we also proposed a A-system that can control the 
secondary triad reality condition. 







(b) An alternative guideline is to obtain as many non-zero eigenvalues as one can. More precisely, this case is 
supposed to have pure imaginary eigenvalues. In such a case, the constraint propagation equations (e.g. 
dtC = ±ikC) behave like the normal wave equations in its original component (e.g. dtC = ±d x C), and 
its stability can be discussed using von Neumann stability analysis. As is well known, stability depends 
on the choice of numerical integration scheme, but it is also certain that we can control (or decrease) the 
amplitude of the constraint terms. 

The advantage of this adjusted system is that we do not need to add variables to the fundamental set, while 
the above first guideline (3a) is the same mechanism which is applied for the A-system. We note that the non-zero 
eigenvalue feature was conjectured in Alcubierre et. at. Jl5[ | in order to show the advantage of the conformally-scaled 
ADM system, but the discussion there is of dynamical equations and not of constraint propagation equations. 

The guideline (3b) is obtained heuristically as we will show in Fig.^| that a system with three zero eigenvalues is 
more stable than one with five. We, however, conjecture that systems with non-zero (or pure-imaginary) eigenvalues 
in their constraint propagation equations have more dissipative features than that of zero-eigenvalue system. This is 
from the von Neumann's stability analysis, evaluating dynamical variables with the finite-differenced quantities. See 
Appendix |b] for more details. 

We remark that adding constraint terms to the dynamical equations is not a new idea. For example, Detweiler |l2| ] 
applied this procedure to the ADM equations and used the finiteness of the norm to obtain a new system. This is also 
one of the standard procedures for constructing a symmetric hyperbolic system (e.g. [p|~p||lC|] ) . We believe, however, 
that the above guidelines yield the essential mechanism for our purpose, to constructing a stable dynamical system. 

In the following subsections and the Appendix |A|, we demonstrate that this adjusted system actually works as 
desired in the Maxwell system and in the Ashtekar system of the Einstein equations, in which above two guidelines 
are applied respectively. 



B. Example 1: Maxwell equations 



1. adjusted system 

We here again consider the Maxwell equations ( |2.9| )-( |2~Tl| ). We start from the adjusted dynamical equations 

d t Ei = ceS k d 3 B k + P t C E + P> l {d J C E ) + Q t C B + ^(fyCn), (3.7) 

d t B t = -ca^djEk + R t C E + rh{djC E ) + S t C B + s^(d 3 C B ), (3.8) 

where P, Q 7 R, S,p, q, r and s are multipliers. These dynamical equations adjust the constraint propagation equations 
as 

d t C E = {diP^CE + P l (d,C E ) + {d % Q l )C B + Q t (d l C B ) 

+ (d tP > l )(d 3 C E ) + p>\d l d ] C E ) + (di^idjCs) + (f'ididjCB), (3.9) 
d t C B = {d t R l )C E + R l (d t C E ) + {d t S l )C B + S'&Cb) 

+(dir ji )(djC E ) + ^{didjCE) + {d^ l ){d 3 C B ) + s^{didjC B ). (3.10) 
This will be expressed using Fourier components by 

» (C E \ _(d l P' + iP% + ikj (d lP > 1 ) - kikjpi* d t Q l + iQ% + ikj (d^) - k t k^ \(C E \ = T (C E \ ,„ . , , 
\C B ) \ d l R t + iR l h + ikj for*) - hkjri* d . S i + + ik . ( 9 . s ii) _ k . kjS ji ) \c B ) \C B ) ' K ' 

Assuming the multipliers are constants or functions of E and B, we can truncate the principal matrix as 



( ) T _ / iP l h - kikjp> % iQ l h - hkjq 31 \ . . 

I ;/,"/,-, - fcifcj-r* iS% - /,•,/,•,>••" ' ' ' 



with eigenvalues 



A ± = p + s±y/p 2 + 4:qr-2ps + s 2 ^ ^ 



where p := iP l ki — kikjp^ 1 , q := iQ % ki — kikjq^ 1 , r := iR % ki — kikj^ 1 , s := iS l ki — kikjS^. 

If we fix q = r = 0, then = p, s. Further if we assume pi" 1 , s^ 1 > 0, and set everything else to zero, then A ± < 0, 
that is we can get the all eigenvalues which have negative real part. That is, our guideline (a) is satisfied. (Conversely, 
if we choose q = r = and p> 1 , s Ji < 0, then A ± > 0.) 
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2. Numerical Demonstration 



We applied the above adjusted system to the same wave propagation problem as in § Q B 3j . For simplicity, we fix 
K = p 1 ^ — s y and set other multipliers equal to zero. In Fig||, we show the L2 norm of constraint violation as a 
function of time, with various k. As was expected, we see better performance for n > (of the system with negative 
real part of constraint propagation equation), while diverging behavior for k < (of the system with positive real 
part of constraint propagation equation). 



o 

E -1.00 



K* -0.1; 



K= + 0.3 



K= + 0.2 



-2.50 I — ' — ' — ' — I — ' ' — ' — I — ' — ' — ' 1 — ' — ' — ' — I — ' ' — ' — 1 

0.0 2.0 4.0 6.0 8.0 10.0 

time 



FIG. 4. Demonstrations of the adjusted system in the Maxwell equation. We perform the same experiments with § [I B 3 
[Figjlj. Constraint violation (L2 norm of Ce) versus time are plotted for various k(= p-'i = s J i). We see k > gives a better 
performance, (i.e. negative real part eigenvalues for the constraint propagation equation), while excessively large positive n 
makes the system divergent again. 



C. Example 2: Einstein equations (Ashtekar equations) 

1. Adjusted system for controlling constraint violations 



We here only consider the adjusted system which controls the departures from the constraint surface. In the 
Appendix, we present an advanced system which controls the violation of the reality condition together with numerical 
demonstration. 

Even if we restrict ourselves to adjusted equations of motion for (J5*, Af) with constraint terms (no their derivatives), 
generally, we could adjust them as 

d t K = -iVj (e cb a NEiED + 2V, (N^E* ) + iA h Q e ab c E l c + X* a C u + Y?C Mj + PfC Gb , (3.14) 
d t Af = -ie ab c NE 3 b FC + N^Ffi + + ANEf + QfC H + Ri ja C Mj + ZfC Gb , (3.15) 

where X l a , Y£ 3 ', Zf b , Pl b ,Qf and are multipliers. However, in order to simplify the discussion, we restrict multipliers 
so as to reproduce the symmetric hyperbolic equations of motion Jl0| , p| , i.e., 

X = Y = Z = 0, Pf = Kl (N' l 5 b a +iNe a bc E t c ), Q? = K 2 (e~ 2 NE?), Rj a = K 3 (i e - 2 Ne ac b E b t Ei). (3.16) 

Here k± = «2 = ^3 — 1 is the case of symmetric hyperbolic equation for (E l a ,A%), while k± — «2 = «3 = is the 
(Ashtekar's original) weakly hyperbolic equation, and other choices of ks let the equation satisfy the level of strongly 
hyperbolic form. 

With these adjusted terms, the constraint propagation equations become 
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d t WC H = (l-2K 3 )(d j WC Mj ), 
d t ^C Mi = (1 - 2 K2 )(d i ^C H ) +iK 3 e mj i(d m ^C Mj ), 
d t {1) C Ga = -2 K3 ^C Ma + t Kl e a bm (d m ^C Gb ). 



(3.17) 
(3.18) 
(3.19) 



against the Minkowskii background. The eigenvalues of the coefficient matrix after the Fourier-transformation are 

(0, ±iKiVk?, ±in 3 V~k?, ±i(2K 2 - 1)(2k 3 - 1)VP) (3.20) 
where k 2 := kik 1 . For example, 



(0 (multiplicity 5), ztivk 2 ) 
(0 (multiplicity 3), ±iV~k? (multiplicity 3)) 



for Kl = K2 = k 3 — 
for Kl = K2 = k 3 = 1 



original system (3-21) 
symmetric hyperbolic system. (3.22) 



That is, our guideline (b) is obtained. 

The above adjustment, (3.14)-(3.16), will not produce negative-real-part eigenvalues, so our guideline (a) cannot be 
applied here. If we adjust the dynamical equation using the spatial derivatives of constraint terms, then it is possible 
to get all negativ e eige n value s like in the Maxwell system (though this is complicated). However, since we found that 
this adjustment, (3.14)-(3.16), gives us an example of controlling the violation of constraint equations for our purpose, 
we only show this simpler version here. 



2. Numerical Demonstration 



As a demonstration, we use here the same model as in §HC 3, that is, gravitational wave propagation in the plane 



symmetric spacetime, with an artificial error in the middle of time evolution. We examine how the adjusted multipliers 
contribute to the system's stability. In Fig.[|, we show the results of this experiment. We plot the violation of the 
constraint equations both Ch and Cmx- An artificial error term was added at t — 6, as a kick of A 2 — ► A 2 (I + error), 
where the error amplitude is +20% as before. We set k = K\ = K2 = ^3 for simplicity. The solid line is the case of 
k = 0, that is the case of "no adjusted" original Ashtekar equation (weakly hyperbolic system). The dotted line is for 
k = 1, equivalent to the symmetric hyperbolic system. We see other line (k = 2.0) shows better performance than 
the symmetric hyperbolic case. 




FIG. 5. Demonstration of the adjusted system in the Ashtekar equation. We plot the violation of the constraint for the 
same model with Fig.^(b). An artificial error term was added at t = 6, in the form of Ay — » A y (l + error), where error is +20% 
as before. Fig. (a) and (b) are L2 norm of the Hamiltonian constraint equation, Ch, and momentum constraint equation, Cmx, 
respectively. The solid line is the case of k = 0, that is the case of "no adjusted" original Ashtekar equation (weakly hyperbolic 
system). The dotted line is for k = 1, equivalent to the symmetric hyperbolic system. We see other line (k = 2.0) shows better 
performance than the symmetric hyperbolic case. 
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IV. DISCUSSION 



With the purpose of searching for an evolution system of the Einstein equations which is robust against perturbative 
errors for the free evolution of the initial data, we studied two "asymptotically constrained" systems. 

First, we examined the previously proposed "A-system" , which introduces artificial flows to constraint surfaces based 
on the symmetric hyperbolic formulation. We showed that this system works as expected for the wave propagation 
problem in the Maxwell system and in Ashtekar's system of general relativity. However, the A-system cannot be 
applied to general dynamical systems in general relativity, since the system requires the base system to be symmetric 
hyperbolic form. 

Alternatively, we proposed a new mechanism to control the stability, which we named the "adjusted system" . This 
is simply obtained by adding constraint terms in the dynamical equations and adjusting the multipliers. We proposed 
two guidelines for specifying multipliers which reduce the numerical errors; that is, non-positive-real-part or pure- 
imaginary eigenvalues of the adjusted constraint propagation equations. This adjusted system was also tested in the 
Maxwell system and in Ashtekar's system. 

As we denoted earlier, the idea of adding constraint terms is not new. However, we think that our guidelines for 
controlling the decay of constraint equations are appropriate for our purposes, and were not suggested before. Up to 
our numerical experiments, our guidelines give us clear indications whether the constraints decay (i.e. stable system) 
or not for perturbative errors, though we also think that this is not a complete explanation for all cases. This feature 
may be explained or proven in different ways, such as finitcness of the norm (of evolution equations or of constraint 
propagation equations), or by another mechanism in future. 

Secondary conclusion is that the symmetric hyperbolic equation is not always the best one for controlling stable 
evolution. As we show in the wave propagation model in the adjusted Ashtekar's equation, our eigenvalue guidelines 
affect more than the system's hyperbolicity. (We found a similar conclusion in We think this result opens a 

new direction to numerical relativists for future treatment of the Einstein equations. 

We are now applying our idea to the standard ADM and conformally scaled ADM system to explain these differences. 
Results will be reported elsewhere jlTj . 
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APPENDIX A: CONTROLLING REALITY CONDITION BY ADJUSTED SYSTEM 

We demonstrate here that our adjusted system in the Ashtekar formulation also works for controlling reality 
conditions. As a model problem, we concern the degener ate point passing problem which we considered previously in 
we review this background briefly, and in § |A 2| we show our numerical demonstrations. 



1. Degenerate point passing problem 

In ]18| , the authors had examined the possibility of dynamical passing of the degenerate point in the spacetime. 
There the authors found that we are able to pass (i.e. continue time evolutions) if we could foliate the time-constant 
hypersurface into complex plane assuming that such a degenerate point exists on the real plane. Such foliations are 
available within Ashtekar's original formulation, since the fundamental variables are complex quantities. The trick is 
to violate the reality condition locally, only in the vicinity of a degenerate point. 

As a model, we construct a metric, ^g, which possesses a degenerate point (det ^g — 0) at the origin t = x = in 
Minkowskii background metric: 

ds 2 = -[1 - (2tacxp(-t 2 - x 2 )) 2 ]dt 2 + 4tacxp(-£ 2 - x 2 )[l - (1 - 2x 2 )exp(-t 2 - x 2 ))dtdx 

+ [1 - (1 - 2a; 2 ) exp(-t 2 - x 2 )fdx 2 + dy 2 + dz 2 . (Al) 
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We consider the time evolution, which initial data is described by a particular time slice t < of (Al), and whose 
time-constant hypersurfaces are foliated by the gauge condition, 



N=l, (N = e- 1 ), 

N x = 2taexp(-i 2 - x 2 )[l - (1 - 2a; 2 ) exp(-i 2 - x 2 )] + iat exp(-b(t 2 + .x 2 )), 
•4g = 0, 



(A2) 

(A3) 
(A4) 



which enables to detour into the complex plane. Our goal is to demonstrate that the time evolution comes back to 
the real plane without any divergence in variables and curvatures. Such a "recovering condition" can be described by 



SsN(t,x)dt = 0, 



?SN l (t, x)dt = 0, Foliation recovering condition 



SSN(t,x) -> 0, SSNtfax) -> 0, 3 E l a E ja /detE(t,x) -> 0, Asymptotic reality condition 
for all four limits x — > ± Ace, t — > t* ± Ai. 



(A5) 
(A6) 



Numerically, this problem becomes an eigenvalue problem, since our boundary conditions, (A5) and (|A6|), specify 



much freedom. To see if the evolution satisfies the criteria or not, we introduced two measures 
F(t fi na i) := max 1 3? (e(t = tr inaU x) - 1)| (asymptotically flat) 

X 

R(t final) := max \5s (e(f = tf ina i,x))\ (asymptotically real) 



(A7) 
(A8) 



and searched the parameters a and b in (A3). 

If we apply our adjusted system to this model, then we expect that the allowed range for the parameters a and b 
becomes more general, since the real-surface-recovering feature is in the flow of the adjusted system's foliation. 



2. Application of the adjusted system 

As was shown in the previous section, for this purpose, we have to foliate our hypersurface in the complex-valued 
region and foliate back to the real-valued surface. That is, we can treat the reality condition, both primary and 
secondary, as a part of the constraint equations. 

For the above degenerate point-passing problem, we nee d to control only the violation of ^sm(E^E^). Therefore, 



similar to the proposal of the adjusted system discussed in § II C , our adjusted dynamical equations can be written as 



d t E l a = -W^NEiEl) + 2V 3 {N^Ej) + iA^JK + X a C H + Y»C Mj + if C G6 + T l Qjfc 3m(^^), (A9) 
d t At = -te ab c N& b F^ + N^F- + V. t A a Q + ANEf + Q a t C H + R?C Mj + ZfC Gb + r tjk 3m(^^), (A10) 

where X*, Y^,Zf b , P* b , Qf, RV ,T l ajk and V a ijk are adjusted multipliers. 

If we simply set X l a = Yj> = Zf = P lb = Qf = Rf j = V a ljk = and T\ jk = -iK5ft ak , (where k is real constant), 
then we obtain the constraint propagation equation 

d t (( ^(EiEi)) = -2K(^(E a El)) + other constraint terms. (All) 

The eigenvalue of the Fourier-transformed RHS is —2k. That is, if we set k > (< 0) then the eigenvalue is negative 
(positive), while k = recovers the original non-adjusted system. 

The results of numerical demonstration are shown in Fig.^. We plot the L2 norm of the violation of the reality 
condition as a function of time, t (this evolution is from t = —5 to 5 ]l8|). Around the time t = the error appears 
due to our "detour" slicing condition, and the original system (k = 0) will not recover the reality surface with the 



choice of a and b in (A3) for this plot. However, for the positive k case, the foliation will be forced to recover the 
reality surface, while for negative k case it will not. 

Therefore this example again supports our guidelines, i.e. negative eigenvalue of constraint propagation equation 
will guarantee the evolution to the constraint surface. 
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time 



FIG. 6. Demonstration of the adjusted system to control the reality condition in the Ashtekar formulation. Reality violation 
(L2 norm of imaginary part of density) versus time are plotted for various adjusted coeffici ent re = —1,0, 1,2. We see re > 
has better performance, (negative real-part eigenvalues of the reality propagation equation, ( All )). 



APPENDIX B: VON NEUMANN ANALYSIS OF CONSTRAINT PROPAGATION EQUATIONS 



Here we show von Neumann's stability analysis for the constraint propagation equations, in order to support our 
guideline (3b) for the adjusted system (§111 A). The von Neumann analysis (see e.g. Jl9[ ) gives us powerful predictions 
for the stability of a finite difference approximation. Briefly, the analysis consists from the Fourier decomposition in 
the spatial directions of the dynamical variables and its one-step time evolution with a particular time integration 
scheme. If we wrote the fundamental variable cf>(x,t), then the criteria for the stability is |Aj| < 1 where A.; are 
the eigenvalues of the amplification matrix G, which is in the expression of the evolution equations in the form of 
<f>{x,t + At) = G<p(x,t). 

In our discussion, the constraint propagation equations are not directly used for numerical integrations, but are used 
as a guideline for the stability. The application of von Neumann analysis, however, is also allowed for the constraint 
propagation equations, as far as substituting the finite derivatives in the analysis using those of the fundamental 
dynamical variables. Here we show the most simplest cases for the adjusted Maxwell system and the adjusted 
Ashtekar system. 

a. Adju sted Maxwell system We start from choosing re := P± = Pi = P % and other multipliers zero in the system 
(3.7) and (3.8). T he Fo urier component of the propagation equation for Ce ( 3.11 ) becomes dtCE — in{k x 



which eigenvalue (3.13) is in{k x + k y 4- k z ). That is, non-zero re gives us a pure-imaginary eigenvalue. By applying von 



Neumann analysis, we obtain the amplitude Gs for FTCS (forward time and center space difference), Brailovskaya 
and (2-iteration) iterative Crank Nicholson schemes as 



|G FT cs| 2 -l + M 2 , 



IG, 



CN2\ 



1 



(rea) 2 + (rea) 4 , 
(rea) 4 /4 + (re<7) 6 /16, 



(Bl) 
(B2) 
(B3) 



respectively, where a = (At / Ax)(sm(k x Ax) + sm(k y Ax) + sin(£; z An)) and we assume 3-dimensional finite grid of 
equal space Ax in all directions. Except for the FTCS scheme, we sec that n on-ze ro |re| (near re = 0) yields |G| < 1. 
The bigger |re| (near re = 0) gives less |G| < 1. The simulation we showed in §111 B is not this case (since we tried to 
show the one which satisfy the guideline (3a)), but we also obtained the numerical results which confirm our conjecture 

here. 

b. Adjusted Ashtekar system Similarly, for the constraint propagation equations (3.17)- ( |3.19| ) 
K2 = K3, we obtain the eigenvalues Xi of the amplification matrix G for the above three schemes, 



with re := rei = 



|2 



l + (recr) 2 , 1 + {(1 - 2re)(re<r)}l 5 



X\ 2 Br = l, 1 - (nay + (rea) 4 , 1 - {(1 - 2re)a}' + {(1 - 2re)a} 4 



|A| 



CN2 



= 1, l-(reCT) 4 /4+(re(j) 6 /16, 1 - {(1 - 2re)cr} 4 /4 + {(1 - 2re)cr} 6 /16, 



(B4) 
(B5) 
(B6) 
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with multiplicity 1, 4 and 2, respectively. Here again we see that non-zero |k| makes the system \G\ < 1 for Brailovskaya 
and 2-iteration Crank-Nicholson schemes. This analysis supports why the guideline (3b) works for our results shown 
in Fig. |. 
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